clear;
close all;

qbit=3;showfig = 0;

Nsub=256;   % no. of subcarriers for rate computation
Ngrid = 1000;
L = 20e-1;
fc = 30e9;
c = 3e8;
BW_des = 1e9;
Ptx = 10;%dBm
Rates = [];
Ptx_lin = 10.^((Ptx-30)/10); % in watts
Lb = c/fc;
spac = Lb/2;
hp = 6.625e-34;
kb = 1.38e-23;
Tk = 290;

factor = 1.6;
x1 = 40e-1*factor;
y1 = 50e-1*ones(size(factor));
z1 = 90e-1*factor;
x2 = -30e-1*factor;
y2 = 50e-1*ones(size(factor));
z2 = 40e-1*factor;
TRX = zeros(6,length(factor));
TRX(1,:) = x1;
TRX(2,:) = y1;
TRX(3,:) = z1;
TRX(4,:) = x2;
TRX(5,:) = y2;
TRX(6,:) = z2;
dist = rms(TRX(1:3,:))*sqrt(3)+rms(TRX(4:6,:))*sqrt(3);

N = round(L/spac);
xcoord=linspace(-L/2,L/2,N);
ycoord=xcoord;
[xx,yy]=meshgrid(xcoord,ycoord);
a = (sqrt((xx-x1).^2+(yy-y1).^2+z1^2)+sqrt((xx-x2).^2+(yy-y2).^2+z2^2))/2;
a_vec = linspace(min(a(:)),max(a(:)),Ngrid)';% column vector
Distmat = a*2;
distgrid = a_vec*2;
Distmat1 = sqrt((xx-x1).^2+(yy-y1).^2+z1^2);
Distmat2 = sqrt((xx-x2).^2+(yy-y2).^2+z2^2);

Dmax = max(Distmat(:));
Dmin = min(Distmat(:));
deltaD = Dmax-Dmin;


showfig = 1;
[ModBFmat,QBFmat,modfn] = Beamforming(qbit,L,TRX,fc,BW_des*1,Ptx,showfig);
freqvec=2*pi*linspace(fc-BW_des/2,fc+BW_des/2,Nsub*2);
antt_delays=(Distmat(:))/c;
AF_product=antt_delays*freqvec;
h_0 = (Distmat1(:)*freqvec/c).*(Distmat2(:)*freqvec/c);

Wb_mch=exp(-1i*AF_product);
SCpowers_new=abs(ModBFmat(:).'*Wb_mch).^2;
SCpowers_old=abs(QBFmat(:).'*Wb_mch).^2;
SCpowers_ideal = zeros(size(freqvec));
SCpowers_ideal(abs(freqvec/2/pi-fc)<=BW_des/2) = 1;
Pow = sum(modfn.^2)*spac^-4*(a_vec(end)-a_vec(1))/Ngrid*c/(freqvec(2)-freqvec(1))*pi;
temp = [Pow sum(SCpowers_old)];
SCpowers_ideal = SCpowers_ideal*Pow/sum(SCpowers_ideal);
SCpowers_ideal = SCpowers_ideal+min([SCpowers_old SCpowers_new])/2;
scale = 1/SCpowers_ideal(round(end/2))*10;
SCpowers_new = SCpowers_new*scale;
SCpowers_old = SCpowers_old*scale;
SCpowers_ideal = SCpowers_ideal*scale;


[ModBFmat,~,~] = Beamforming(qbit,L,TRX,fc,BW_des*1.2,Ptx,showfig);
SCpowers_new_2 = abs(ModBFmat(:).'*Wb_mch).^2;
SCpowers_new_2 = SCpowers_new_2*scale;

[ModBFmat,~,~] = Beamforming(qbit,L,TRX,fc,BW_des*1.3,Ptx,showfig);
SCpowers_new_3 = abs(ModBFmat(:).'*Wb_mch).^2;
SCpowers_new_3 = SCpowers_new_3*scale;

        Sinc = sinc((freqvec/2/pi-fc)*deltaD/c);
        SCpowers_sinc = N^4*Sinc.^2*scale;
% 
% figure(1);
% hold on
% plot(freqvec/2/pi/1e9,10*log10(SCpowers_old)-max(10*log10(SCpowers_old)),'--','LineWidth',1.5,'MarkerSize',1.5)
% plot(freqvec/2/pi/1e9,10*log10(SCpowers_sinc)-max(10*log10(SCpowers_old)),'*-','LineWidth',1.5,'MarkerSize',1.5)
% hold off
% xlabel('Frequency (GHz)')
% ylabel('Beamforming gain (dB)')
% % title(sprintf('Channel Gain\n BW = %.1f GHz',BW_des/1e9))
% ylim([-50 10])
% legend('Narrow band beamforming','Approximated gain', 'FontSize', 16)
% hold off
% % set(get(gca, 'Legend'), 'Box', 'off'); 
% set(gca, 'FontSize', 16);
% set(gcf, 'Position', [100 100 600 500]); % 调整当前图窗的大小
% set(0, 'DefaultAxesFontName', 'Times New Roman');
% saveas(gcf,'fig/Fig_beamforming_performance_sinc.fig')
close all
figure(2)
hold on
plot(freqvec/2/pi/1e9,(SCpowers_old)/max((SCpowers_old)),'b-','LineWidth',2,'MarkerSize',1.5)
plot(freqvec/2/pi/1e9,(SCpowers_sinc)/max((SCpowers_sinc)),'r--','LineWidth',2,'MarkerSize',1.5)
hold off
grid on
xlabel('Frequency (GHz)')
ylabel('Normalized beamforming gain')
% title(sprintf('Channel Gain\n BW = %.1f GHz',BW_des/1e9))
ylim([0 1.2])
legend('Narrow band beamforming','Approximated gain', 'FontSize', 16,'Location','northwest')
hold off
set(get(gca, 'Legend'), 'Box', 'off'); 
set(gca, 'FontSize', 16);
set(gcf, 'Position', [100 100 600 500]); % 调整当前图窗的大小
set(0, 'DefaultAxesFontName', 'Times New Roman');
saveas(gcf,'fig/Fig_beamforming_performance_sinc1.png')

